Entanglement perturbation theory for the elementary excitation in one dimension 
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The entanglement perturbation theory is developed to calculate the excitation spectrum in one di- 
mension. Applied to the spin-^ antiferromagnetic Heisenberg model, it reproduces the des Cloiseaux- 
Pearson Bethe ansatz result. As for spin-1, the spin-triplet magnon spectrum has been determined 
for the first time for the entire Brillouin zone, including the Haldane gap at k — n. 
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The importance of elementary excitations in condensed 
matter systems may be best understood in the super- 
fluid 4 He. The Tisza two-fluid model with the experi- 
mentally found phonon-roton spectrum explains funda- 
mental properties of the superfluid 4 He Feynman's 
effort then to explain the roton spectrum is well known 
Q. From the theorem of Bloch-Floquet, the elementary 
excitation with momentum k for a translationally invari- 
ant Hamiltonian H is written as 



/V 



* fc = 5> ttI Oj|p) 



(1) 



i=i 



where \g) is the ground state and the summation over I 
extends over the entire lattice sites. The Oi is a local 
cluster operator to be determined for a given Hamilto- 
nian. 

In spite of a simplicity and validity of the expression 
((T|) , not much progress has been made along this line since 
the days of Feynman. The Heisenberg antiferromagnet 
(HA) described by the Hamiltonian 



H 



J Si ■ Sj-| 



(2) 



is probably the best studied system concerning the exci- 
tation spectrum. In particular, Haldane conjectured in 
1983 that the half-odd integer and integer spins might 
behave essentially differently [3j, which together with a 
field theoretic prediction of Affleck [4j for a logarithmic 
correction to the power-law behavior in the spin-spin cor- 
relation function in the spin-i case, triggered an inten- 
sive study of HA ranging from the exact diagonalization 
d, @] and Monte Carlo 7,J1 to DMRG (density matrix 
renormalization group) 9h11|. These studies al ong with 
the Bethe ansatz solution for the spin-i case [12[ lead 
to a confirmation of the both claims. Concerning the 
elementary excitation for the entire Brillouin zone, how- 
ever, Takahashi's two attempts following the Feynman 
variational method for 4 He and a projector- Monte Carlo 

, [l4| . And none of the 
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method were the only studies 
previous studies gave a serious consideration to the ex- 
pression (J]). 



In this Letter, we analyze (Q} exactly for the HA ([2]) 
with periodic boundary conditions by the recently devel- 
oped entanglement perturbation theory (EPT). EPT is a 
novel many-body method which takes into account corre- 
lations systematically. Its mathematical implementation 
is singular value decomposition (SVD), intuitively divide 
and conquer. EPT has addressed so far classical statis- 
tical mechanics [Hj], ID quantum ground states [3] and 
2D quantum ground states [TtJ ■ We here address the el- 
ementary excitation in one dimension. By EPT, we are 
not only free from a negative sign problem which is inher- 
ent to MC for quantum spins and fermions, but can also 
handle an order of magnitude larger systems than MC 
and DMRG. The key of the success lies in our ability of 
calculating the ground state \g) precisely and most im- 
portantly in an un-renormalized form. We examine the 
cluster excitation operator Oi systematically. We have 
found that the size of the magnon is 4 lattices long at 
largest for both spin- 1 and spin-1. 

We solve the problem ([T]) in two steps. First, we find 
the ground state. We can do this in two ways. Ei- 
ther to consider the density matrix eigenvalue problem 
e^ ff |a) = e~ pE <>\g) with /3 -> as done in |16] (EPT- 
gl). There, the density matrix is expressed in a matrix, 
tensor product form, reducing the problem to that of the 
partition function calculation in the 2D, 3D Ising model 
[15| . This method is particularly suited to study the in- 
finite system because then we only need to consider the 
largest eigenvalue. Or to directly minimize the energy 
(g\H\g)/(g\g) as formulated in [l7| (EPT-g2). The sec- 
ond method is free from the parameter /3, but it needs 
to handle a number of eigenvalues no matter how large 
the system size. Either way, the starting wave function 
is that of SVD-ed one [H-[l3], 

lfl) = ---Ci / j(«i)C| 7 (*a)^(«3)--- (3) 

on the basis • ■ • |si) ® [si) (g) \sg) ® ■ • • where \si) denotes 
spin states on the site i, and the local bi-partite wave 
functions C, 1 ' 2 reflect the local antiferromagnetic interac- 
tion. 

Let us consider the variational problem, EPT-g2, 
\9\H\g) 



SE a 



(g\g) 



o 



(4) 
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size 


EPT 


BA 


16 


-0.4463935 


-0.4463935 


64 


-0.4433459 


-0.4433485 


256 


-0.4431555 


-0.4431597 



TABLE I: Comparison of the ground state energies between 
BA and EPT with the entanglement p — 15 for 16 spins, 22 
for 64 and 36 for 256 
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FIG. 1: (a) Matrix A. The vertical lines denote entanglement 
while the horizontal ones the basis spin states. The upper 2 
external verticle lines constitute the right index of A while the 
lower 2 the left index. The inner lines should be contracted, 
(b) Matrix A 



with respect to C 1 ' 2 . First consider S(g\g). For the peri- 
odic ring of N bi-lattice units, we have 



(g\g) = T r A 



N 



(5) 



with the matrix A shown in Fig. la. The variation with 
respect to C 1 ' 2 is, due to the translational symmetry, N 
times the variation of £ 1:2 at a particular bi-lattice unit. 
The idea is to use an initial trial for £ 1,2 and evaluate 
the trace over the TV — 1 bi-lattice units, which can be 
done as follows. Writing the left and right eigenvector 
matrixes as L and R and the diagonal eigenvalue matrix 
as a, the matrix A is written as 



A = RaL 



(6) 



where L is the transpose of L, with the property L-R = 1. 
We can thus write 



(7) 



Now note that the last form ([7|) contains C, 1,2 quadrat- 
ically, thus written as either C Ml {QC l or <f 2 M 2 (C)C 2 
with 



Sls,rs(?i 1 L i (lt,rt)R t (lu,ru) 



(8) 



where the indexes are put together as I = (lu — 1) • p ■ 
n s + (Id ~ 1) • n s + Is and r — (ru — 1) • p ■ n s + (rd — 
1) • n s + rs where p denotes the entanglement size and 
n s the local spin degrees of freedom, and summations 
are implied for the repeated indexes. And likewise for 
M 2 ((). As for the numerator (g\H\g), because of the 
translational symmetry, the variation with respect to C, 1,2 
is again TV times the variation of £ 1,2 at a particular bi- 
lattice unit. Since the Hamiltonian is a sum of local 
nearest neighbor interactions, a typical term in (g\H\g) 



distance 


EPT 


BA 


1 


-0.1477187 


-0.1477157 


2 


0.0606790 


0.0606798 


3 


-0.0502424 


-0.0502486 


4 


0.0346217 


0.0346528 


5 


-0.0308335 


-0.0308904 


6 


0.0243932 


0.0244467 


7 


-0.0224726 


-0.0224982 



TABLE II: Comparison of G z over the first 7 sites between 
BA and EPT with p = 36 for 256 spins 



is LjARiVj 2 LiARj(jf 1 where A is given by Fig. lb. A 
summation is carried out over the terms containing £ 1,2 
except a quadratic term to be variated, and (g\H\g) can 
be written as either (f^HOC 1 or C 2 A f2 (C)C 2 - Thus the 
variational problem Q leads to nonlinear, generalized 
eigenvalue problems 



M\QC =E g N i {0C 



i= 1,2 



(9) 



We have solved © iteratively for both spin-i and spin-1 
for the system size up to 1024 and the entanglement size 
up to 38. A thorough discussion of the EPT result for the 
ground state properties of the xxz model with the Ising 
anisotropy coupling A = ~ oo will be given in a future 
publication [18(. We here concentrate on the isotropic 
xxx case. A convergence in the ground state energy per 
site occurs typically around p=20. 

For spin- i, the exact result is given by Bethe ansatz 
(BA) [19(. The ground state energies per site for 16, 64 
and 256 spins by EPT and BA are compared in Table I. 
The spin-spin correlation function G z (l) = (S z (0)S z (l)) 
is also known by BA up to lattice separations 7. Com- 
parison of BA and EPT for 256 spins is given in Table 
II. As for a longer distance I > 1, an asymptotic for- 
mula is given by field theory Extensive finite size 
analysis with exact diagonalization, Monte Carlo and 
DMRG [1, 0, [ll[ have been done concerning the field 
theoretic prediction of the logarithmic correction of the 
form oc y/ln(l)/l at large distance. The lattice size at 
which the asymptotic formula is realized is estimated to 
be several thousands. Fig. 2 shows the EPT result of the 
correlation function for the spin-^ with 256 spins which 
is a lot longer than previously studied. Note that the 
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FIG. 2: Log-log plot of the spin-spin correlation function for 
256 spins for spin-i. The upper line is the converged result 
for entanglement 32, while the lower line is the asymptotic 
formula from field theory. 



converged and hence exact result at p — 32 is well fitted 
by a straight line, except near the central part at I = 128 
reflecting the periodic boundary condition, 



G z (l) oc I 



-0.92 



(10) 



While the greater than -1 exponent reflects the logarith- 
mic correction, a large discrepancy between the EPT and 
the field theoretic asymptotic form, FT in Fig. 2, for a 
distance I ~ 100 is consistent with the argument that 
the FT asymptotic form is correct for distances as large 
as several thousands. It should be noted that EPT can 
handle any system size, but for a larger system, a larger 
entanglement is necessary to get a good ground state. 
An interesting question then is, to calculate the corre- 
lation functions correctly for thousands-long distances, 
how large the entanglement p should be? In fact the is- 
sue is closely related to another field theoretic prediction 
[3| that the correlation functions at large distance should 
have a singular dependence on the Ising anisotropy pa- 
rameter A, namely G z — G x at A = 1 but G z — 4G X 
when A approaches 1 from below. Following the success- 
ful EPT analysis of such symmetry breaking in the 2D, 
3D Ising models 15| , a calculation is currently underway 
and will be reported elsewhere. 

As for spin-1, EPT agrees with DMRG, e.g. the ground 
state energy per site for N = 48 is -1.401482 (EPT with 
p = 24) vs -1.401484 (DMRG 9]), and the spin-spin cor- 
relation function shows an exponential decay as expected 
for a gapped system. 



An important note on the ground state algorithms 
EPT-gl,2 is that they can also calculate some excited 
states. For example, we have applied EPT-g2 to spin-1 
to get the energy gap at k = ir correctly, 0.4124 (EPT 
with p = 20) vs 0.4123 (DMRG @) for 48 spins. 

We now come to the second step, the implementation 
of the variational program for the elementary excitation 
(EPT-e), 



SE k = S^ H l^=0 
Using ([T]), we first rewrite E k as 



l=i 



(11) 



(12) 



where [• • • ] means a commutator and we have used the 
fact that H\g) = E g \g), and therefore the accurate 
ground state wave function is a crucial ingredient in this 
method. Note that N is the number of bi-lattice units 
and 2N is the total number of spins. More importantly, 
we need the un-renormalized ground state which may not 
be easy to obtain by DMRG without an efficient restora- 
tion of the un-renormalized ground state at the end of the 
calculation, particularly for large systems such as 512 and 
1024 spins as studied by EPT. 

To carry out the variation Se(k) — with respect to 
the cluster operator Oi , the simplest case is the one which 
acts only on one site and, say, for the spin-triplet excita- 
tion, then it is uniquely S + = S x + iS y for spin-i. In the 
spin-1 case, there are only two such operators, [i\ = S + 
and H2 = S + S z . The cluster operator can then be writ- 
ten as a linear combination Oi — ci//i + C2/J.2 and the 
variation is with respect to the vector x — {01,02), lead- 
ing to a generalized eigenvalue problem 



Tx = e(k)Ux 



(13) 



where T and U are 2x2 matrixes. The n-cluster op- 
erator Oi is generally written as a linear combination of 
operator products of n n s x n s local operators, where 
n s = 2 for spin-i and 3 for spin-1. With the increase 
of the cluster size of Oj, the matrix size of T and U in- 
creases like 1,4,15 and 56 (spin-i) and 2,16,126 and 1016 
(spin-1) for the spin-triplet excitation. For example, for 
n = 2 and spin-^, Oi is a linear combination of the 4 
local excitation operators, S + <£> 1, S + ® S z , 1 (g> S + and 
S z (g> S + . The calculation of the matrixes T and U arc 
essentially the same as in the ground state, although due 
to the cluster nature of the operator Oi and the pres- 
ence of the commutator [H, Oj] , we have a little lengthy 
algebraic procedure, details to be presented elsewhere. 

Fig. 3 shows the spin-triplet excitation spectrum for 
spin- i with the chain size 512 and the cluster size n up 
to 4. The EPT calculation almost converged at n = 4, 
and gives an agreement of 1% precision with BA (l2^ . 
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FIG. 3: Spin triplet excitation spectrum for spin-i for 512 
spins. The ground state used is from p=20. From top to 
bottom, the cluster size is 1 to 4 (red online). The thick line 
(blue online) is the Bethe ansatz result [l2| . 



FIG. 4: Spin triplet excitation spectrum for spin-1 for 512 
spins. The ground state used is from p=24. From top to 
bottom, the cluster size is 1 to 3 (red online). 



Fig. 4 shows the same for spin-1 up to the cluster size 
n = 3 where the calculation almost converged. The Hal- 
dane gap at k = ir is found by EPT to be 0.414 agreeing 
with the previous results [1, [l(J EH ■ As for the region 
^ k/n ^ 1/2, the spin-triplet spectrum is believed to be 
embedded in a continuum spectrum of a pair of magnons 
with the total spin-z component to be [2(J. The exact 
diagonalization for N = 14 indeed tells us that the lowest 
excitation at k = is spin singlet, presumably a (— n, 7r) 
pair of spin-triplet magnons from k = ±7r [g, |l3[ . We 



believe that the spin-triplet magnon spectrum for spin-1 
for the entire Brillouin zone has been determined for the 
first time by EPT. Moreover EPT can handle not only 
the lowest but the entire spin-triplet sectors, and can be 
repeated for other excitations such as spin-singlet. It 
simply amounts to calculating not just the minimum but 
all the eigenvalues of (fT3")) . Such entire excitation spectra 
have been known only for the isotropic spin-^ case by 
Bethe ansatz, see Fig. 4, 5 in (l9j . 

In conclusion, we have developed EPT for the elemen- 
tary excitation (JIJ, where the un-renormalized ground 
state \g) plays a central role. A challenge to EPT is 
to calculate the correlation functions for a distance long 
enough to compare with field theory. While the applica- 
tions of EPT-e to ID fermions and bosons are straight- 
forward, we need yet to see how it works in two dimen- 
sions. Finally, the successful calculation of the excitation 
spectrum indicates that EPT can handle various nano- 
structures embedded in correlated host materials, open- 
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ing a possible new look at the Kondo effect [21 1. 
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